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ANALYSIS  OF  AN  HP- NON-CONFORMING 
DISCONTINUOUS  GALERKIN  SPECTRAL  ELEMENT  METHOD 
FOR  WAVE  PROPAGATION* 

TAN  BUI-THANH  t  AND  OMAR  GHATTAS  tt§ 

Abstract.  We  analyze  the  consistency,  stability,  and  convergence  of  an  hp  discontinuous 
Galerkin  spectral  element  method.  The  analysis  is  carried  out  simultaneously  for  acoustic,  elastic, 
coupled  elastic-acoustic,  and  electromagnetic  wave  propagation.  Our  analytical  results  are  developed 
for  both  conforming  and  non-conforming  approximations  on  hexahedral  meshes  using  either  exact 
integration  with  Legendre- Gauss  quadrature  or  inexact  integration  with  Legendre- Gauss- Lobat to 
quadrature.  A  mortar-based  non-conforming  approximation  is  developed  to  treat  both  h  and  p 
non-conforming  meshes  simultaneously.  The  mortar  approach  is  constructed  in  such  a  way  that 
consistency,  stability,  and  convergence  analyses  for  non-conforming  approximations  follows  the  con¬ 
forming  counterparts  with  minimal  modifications.  In  particular,  sharp  hp-convergence  results  are 
proved  for  non-conforming  approximations  for  time  dependent  wave  propagation  problems  using 
inexact  quadrature. 

Key  words,  discontinuous  Galerkin  method;  spectral  element  method;  linear  hyperbolic  equa¬ 
tions;  acoustic,  elastic,  and  electromagnetic  wave  propagation;  Riemann  flux;  non-conforming  meshes; 
Legendre-Gauss;  Legendre-Gauss-Lobatto;  consistency,  stability,  convergence. 

AMS  subject  classifications.  65N35,  65N12,  65N15 


1.  Introduction.  The  discontinuous  Galerkin  (DG)  method  was  originally  de¬ 
veloped  by  Reed  and  Hill  [21]  for  the  neutron  transport  equation,  but  has  been  ex¬ 
tended  to  other  problems  governed  by  partial  differential  equations  (PDEs)  [9].  In 
particular,  it  has  emerged  as  a  particularly  attractive  method  for  hyperbolic  PDEs 
[8,  10].  Among  its  many  advantages  over  classical  finite  volume  and  finite  element 
methods,  it  has  the  ability  to  treat  solutions  with  large  gradients  including  shocks,  it 
provides  the  flexibility  to  deal  with  complex  geometries,  and  it  is  highly  parallclizable 
due  to  its  compact  stencil.  Perhaps  its  most  important  advantage,  however,  is  its 
ability  to  support  hp  adaptivity  (where  h  refers  to  mesh  size,  and  p  to  local  poly¬ 
nomial  order)  in  a  natural  manner  [5].  Together,  these  advantages  make  DG  a  very 
desirable  method  for  parallel  solution  of  large-scale  hyperbolic  problems  on  adapted 
meshes  (e.g.,  [27]). 

The  question  of  how  to  treat  non-conforming  interfaces  between  elements  due 
to  both  local  p-refinement  and  local  /i-refinement  can  be  addressed  in  several  ways. 
The  non-conforming  approach  of  Kopriva  [14,  18]  replaces  a  non-conforming  face  by 
mortars  that  connect  pairs  of  contributing  elements.  The  actual  computations  are 
performed  on  the  mortars  instead  of  the  non-conforming  faces,  and  the  results  are 
then  projected  onto  the  contributing  element  faces.  Since  the  mortar  approach  main¬ 
tains  the  compact  stencil  of  the  DG  method,  it  makes  adaptivity  highly  parallclizable 
[27].  Moreover,  discrete  stability  and  optimal  convergence  rates  have  been  numerically 
observed  in  practice  [18,  27].  Since  then,  there  have  been  no  attempts  to  theoretically 
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study  the  stability  and  convergence  of  Kopriva’s  mortar-based  non-conforming  ap¬ 
proximation  in  the  context  of  discontinuous  Galerkin  spectral  element  methods.  This 
is  the  subject  of  the  present  article. 

The  celebrated  Lax-Richtmyer  equivalence  theorem  [20]  has  far  reaching  conse¬ 
quences  in  numerical  analysis,  so  much  so  that  it  is  sometimes  called  the  fundamental 
theorem  of  numerical  analysis.  The  theorem  says  that  for  well-posed  linear  differen¬ 
tial  problems,  consistency  and  stability  of  a  difference  method  imply  its  convergence. 
However,  typical  applications  of  the  Lax-Richtmyer  theorem  provide  an  upper  bound 
on  error  that  grows  exponentially  in  time.  In  practice,  it  is  observed  that  the  error 
grows  at  a  much  lower  rate.  Rather  than  rely  on  the  equivalence  theorem,  it  has  been 
shown  that  a  direct  proof  of  convergence,  with  an  error  bound  that  grows  at  most 
linearly  in  time,  is  possible  for  a  class  of  discontinuous  Galerkin  methods  [12,  13]. 

In  this  paper,  we  study  theoretically  the  consistency,  stability,  and  convergence 
of  a  discontinuous  Galerkin  spectral  element  method  (DGSEM)  using  such  a  direct 
proof.  In  particular,  we  present  a  stability  proof  using  an  energy  method  for  the 
DGSEM  with  the  mortar-based  non-conforming  approximation  of  Kopriva  [14,  18]. 
We  expect  that  the  results  of  our  study  can  be  applied  to  a  large  class  of  linear 
conservation  laws.  However,  for  concreteness,  the  proof  is  simultaneously  carried  out 
for  elastic,  acoustic,  coupled  elastic-acoustic,  and  electromagnetic  wave  equations,  as 
exemplary  conservation  laws  governed  by  linear  hyperbolic  PDEs.  Instead  of  using 
an  exact  numerical  quadrature  for  the  mortars  as  in  [18,  27],  for  which  we  are  not 
able  to  prove  stability  in  the  interesting  case  of  Legendre-Gauss-Lobatto  quadrature, 
we  propose  to  employ  a  particular  kind  of  quadrature  rule  that  not  only  facilitates 
the  stability  proof,  but  is  also  cheaper. 

We  are  able  to  prove  stability  and  convergence  for  both  exact  numerical  quadra¬ 
ture  using  the  Legendre-Gauss  rule  and  inexact  numerical  quadrature  using  the  Legendre- 
Gauss-Lobatto  rule.  For  the  inexact  numerical  integration,  we  use  the  tensor  product 
quadrature  rule,  since  it  allows  us  to  perform  discrete  integration  by  parts  [24,  17], 
which  in  turns  paves  the  way  for  the  stability  and  convergence  proofs. 

The  Riemann  numerical  flux  is  our  main  ingredient  in  proving  stability  and  con¬ 
vergence.  Since  the  PDEs  in  this  paper  are  linear,  the  exact  Riemann  flux  can  be 
derived  [25],  and  therefore  we  use  it  in  our  proofs.  As  will  be  shown,  it  is  the  dis¬ 
sipative  nature  of  the  Riemann  flux  that  makes  the  DGSEM  stable.  Therefore,  we 
speculate  that  our  results  also  hold  for  other  dissipative  fluxes  such  as  the  stabilized 
Lax-Friedrichs  numerical  flux  [22]. 

We  mostly  restrict  ourselves  to  the  case  of  affine  hexaliedral  meshes,  for  which 
details  of  the  derivations  and  proofs  are  presented.  In  order  to  make  the  proofs 
concrete,  three-dimensional  coupled  elastic-acoustic  and  electromagnetic  waves  are 
used  as  examples.  Of  course,  all  of  the  results  hold  for  two-dimensional  problems  as 
well. 

This  article  is  organized  as  follows.  Section  §2  briefly  describes  a  weak  setting  for 
general  linear  conservation  laws.  Section  §3  presents  an  hp  DGSEM  for  both  elastic- 
acoustic  and  electromagnetic  waves.  We  then  prove  stability  for  conforming  meshes 
in  §4.  The  detailed  description  of  our  mortar-based  non-conforming  approximation  is 
given  in  §5,  which  is  followed  by  the  non-conforming  stability  proof  in  §6.  The  direct 
proof  of  convergence  is  carried  out  in  §7,  and  Section  §8  offers  conclusions. 

2.  General  setting  for  linear  hyperbolic  conservation  laws.  We  are  in¬ 
terested  in  linear  wave  equations  governed  by  linear  hyperbolic  conservation  laws.  In 
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the  strong  form,  a  general  equation  is  given 


as 


Qvw  +  Vc  •  (3q)  —  0,  c\GV,xGT>, 
at 

with  V  as  the  solution  space,  to  be  specified  later,  over  the  domain  of  interest  V, 
and  with  appropriate  initial  and  boundary  conditions.  The  subscript  x  denotes  the 
^-coordinate  system  in  which  the  divergence  operator  acts.  Next,  multiplying  by  the 
test  function  p,  the  corresponding  weak  formulation  is  obtained  as 


[Q * 

b  dt 


p  dx 


tv 


(3q)  •  p  dx  =  /  0  •  p  dx, 

Jv 


where  denotes  the  Euclidean  inner  product. 

We  next  partition  the  domain  T>  into  Ne\  non-overlapping  hexahedral  elements 
such  that  V  =  U^i  De>  and  integrate  the  weak  formulation  by  parts  twice  to  obtain 

Y^fo  Qe%-^dx  +  JD  V»-(3qe)-p edx 

+  [  n  •  [(5qe)*  -  3qe]  •  Pe  dx  =  V]  /  ge  •  pe  dx,  (2.1) 

Jd  D'  e  J  D« 


where  a  consistent  numerical  flux  (5qe)*  has  been  introduced  to  couple  solutions 
of  neighboring  elements,  and  (-)e  denotes  the  restriction  on  the  e-th  element  of  the 
corresponding  quantity. 

Equation  (2.1)  is  known  as  the  strong  form  in  the  context  of  nodal  discontinuous 
Galerkin  methods  [13].  For  the  DGSEM  described  in  this  paper,  the  strong  form 
(integrating  the  flux  terms  by  parts  twice)  and  the  weak  form  (integrating  the  flux 
terms  by  parts  once)  are  equivalent  [24,  17],  and  hence  all  the  results  in  the  paper- 
hold  for  the  weak  form  as  well. 


3.  A  discontinuous  Galerkin  spectral  element  method.  In  this  section,  we 
briefly  describe  an  /ip-discontinuous  Galerkin  spectral  element  method.  We  approxi¬ 
mate  each  element  De  by  polynomials,  again  denoted  by  De,  such  that  each  element 
De  is  mapped  to  the  reference  hexahedron  D  =  [— 1,  l]3  by  a  C1-diffeonrorphism  Xe, 
and  V  «  VNel  =  U^i  De.  Equation  (2.1)  can  be  written  in  terms  of  D  as 


/9D 


-frf 


■  pe  dr  =  /.  jese  ■  pe  dr’ 


(3.1) 


where  r  =  {r\,r2,T^)  £  D  represents  the  reference  coordinates  and  Je  is  the  Jacobian 
of  the  transformation.  The  outward  normal  on  the  boundary  of  the  master  element 
D  is  denoted  by  n,  and  the  contravariant  flux  [15]  is  defined  as 


$i  =  Jeai-$,  i  =  1,2,3, 


with  a1  as  the  contravariant  basis  vectors. 

We  now  describe  the  approximation  spaces  for  wave  propagation  in  elastic,  acous¬ 
tic,  and  coupled  elastic-acoustic  media  using  the  strain-velocity  formulation,  and  for 
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Maxwell’s  equations  for  electromagnetic  wave  propagation.  Equation  (3.1)  can  be 
specialized  to  the  elastic-acoustic  wave  propagation  case  by  the  following  definitions, 


q  = 


G  V, 


Q 


G  V, 


with  I  denoting  the  fourth-order  identity  tensor,  0  the  zero  tensors  of  appropriate 
sizes,  I  the  second-order  identity  tensor,  E  the  strain  tensor,  v  the  velocity  vector,  / 
the  external  volumetric  forces,  and  p  the  density. 

The  action  of  the  flux  operator  5  on  the  strain-velocity  unknowns  q  can  be  shown 
to  be  [27] 


(ffq)i 


(~\  [y®er  +  e*  (E>  v)\ 

V  -(CJ3)ei  )GV 


for  i  =  l,2, 3. 


For  isotropic  linear  elasticity,  the  strain  tensor  E  and  the  Cauchy  stress  tensor  S  are 
related  by  the  fourth-order  constitutive  tensor  C, 


S  =  CE  =  \tr(E)I  +  2pE, 


where  A  and  p  are  the  two  Lame  constants  characterizing  the  isotropic  constitutive 
relationship.  The  longitudinal  wave  speed  cp  and  shear  wave  speed  cs  are  defined  in 
terms  of  the  Lame  constants  and  density  by 


Cp  — 


'A +  2  [i 


with  cs  =  0  in  acoustic  regions  by  virtue  of  p  =  0. 

As  in  [27],  we  choose  the  solution  space  to  be  V  =  VfyJ[  ®  V 3,  where  V  denotes 
a  space  of  sufficiently  smooth  functions  defined  on  T>  so  that  (2.1)  makes  sense.  The 
discontinuous  approximation  to  V  is  given  by 

Vd  ■=  {qd  G  L\VN°')  :  qd\Dc  o  Xe  G  QNc(D)}, 


where  Qjve  is  the  tensor  product  of  one-dimensional  polynomials  of  degree  at  most  Ne 
on  the  reference  element.  It  should  be  pointed  out  that  the  polynomial  orders  need 
not  be  the  same  for  all  directions.  Nevertheless,  we  use  the  same  order  for  clarity  of 
the  exposition.  The  numerical  solution  q^  €  ®  Vd  restricted  on  each  element 

De  is  specified  as 


qd|Dc  o  Xe  G  Vd  =  Q^fsym  ®  Xe  :  D  DC 

Before  introducing  the  Riemann  flux,  let  us  recall  the  following  standard  DG  notation 
for  quantities  associated  with  element  interfaces: 

[qj  =  q+.n+  +  q-.n-  [q]  =  q^  -  q+,  M  \  Z~  , 

where  the  positive  and  negative  signs  indicate  element  interior  and  exterior,  respec¬ 
tively. 

For  linear  conservation  laws  one  can  solve  the  Riemann  problem  exactly  by  various 
methods  [25].  Using  the  Rankine-Hugoniot  approach,  Wilcox  et  al.  [27]  show  that 
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the  exact  Riemann  flux  for  the  strain  equation  is  given  by 

n  ■  [(£q)js  -  (5q )b]  =  (fco n  •  [CJ5]  +  k0p+c+  [u])  n  ®  n 

—  k\  sym  (n®  (n  x  (n  x  [Ci3])) 

—  k\  p+c+  sym  (n®  (n  x  (n  x  [i>])) , 


and  for  the  velocity  equation  by 

n  '  [(3q)*  -  (5"q)„]  =  (k0n  •  [C£7]  +  k0p+c+  [«])  p~c~n 

—  k\p~c~n  x  (n  x  [C.E]) 

—  kip+cf  p~ c~n  x  (n  x  [u]), 

with  ko  =  ( P~c~  +  p+Cp  )  1 ,  fci  =  (p~ c~  +  p+c+ )  1  if  p~  ^  0,  and  k\  =  0  if  p~  =  0. 
Here,  we  will  consider  only  traction  boundary  conditions  Sn  =  tbc,  where  tbc  is 
the  prescribed  traction.  The  traction  condition  is  enforced  by  the  following  mirror 
principle, 


M  =  M  =  o, 


and  [5]  =  -2  (tbc  -  S  n)  , 


which  applies  to  both  elastic  and  acoustic  media. 

We  next  specialize  (3.1)  to  the  case  of  Maxwell’s  equations.  In  this  case, 


q  = 


G  V,  Q  = 


el  0 
0  pi 


fl  = 


G  V, 


where  E  denotes  the  electric  field,  H  the  magnetic  field,  p  the  permeability,  and  e 
the  permittivity. 

The  action  of  the  flux  operator  ^  on  the  electromagnetic  field  q  is  defined  by 


(Sq)i  = 


— e*  x  H 

e,  x  E 


GV,  for  i  =  1,2,3, 


where  V  =  V3  and  qd|De  o  Xe  G  VJ  =  Q3N  .  The  exact  Riemann  numerical  flux  for 
electric  equation  can  be  shown  to  be  [13] 


[(3q).E  (3q)f;  — 


1 


2  {{%}} 


n  x  (Z+  [H]  —  n  x  [E])  , 


and  for  the  magnetic  equation, 

n  ■  [(ffq)«  -  (ffq)«]  =  -  x  (y+  lE]  +n  x  W)  , 

where 


Similar  to  the  elastic-acoustic  coupling  case,  we  use  the  mirror  principle  to  enforce 
a  perfect  electric  conductor  (PEC)  boundary  condition  by 

Z~  =  Z+,  Y~  =Y+,  nx[H}  =  0,  n  x  [E]  =  2n  x  £T, 
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and  a  perfect  magnetic  conductor  (PMC)  boundary  condition  by 

Z~=Z+,  Y~=Y+,  nx[E\=  0,  n  x  [H]  =  2n  x  H  . 

For  dielectric  boundary  conditions  we  use 

n  x  E~  =  n  x  E+ ,  n  x  H  =  n  x  H  '  . 

In  order  to  unify  the  treatment  for  elastic,  acoustic,  coupled  elastic-acoustic,  and 
electromagnetic  waves,  we  define  a  generic  polynomial  space  Vn  whose  meaning  will 
be  clear  in  each  context.  For  example,  if  we  write  qe  €  Vn,  this  identifies  Vn  =  VJ. 

The  tensor  product  basis  for  Qn  is  built  upon  the  following  one-dimensional 
Lagrange  basis 


m 


TT 

k=0,l,...,N  ^  ^ k 

k^l 


where  the  TVth-degree  Legendre-Gauss-Lobatto  (LGL)  points,  or  ATth-degree  Legendre- 
Gauss  points  (LG),  {£/}  on  [—1,1]  for  l  =  0, . . . ,  N,  are  chosen  as  both  the  interpo¬ 
lation  and  quadrature  points.  This  is  also  known  as  the  collocation  approach.  The 
Lagrange  interpolant  of  a  function  /(r)  on  the  reference  element  D  is  defined  through 
the  interpolation  operator  In  as 


N 

XN  (/)  =  E  flmn^(ri)im(r2)L(r3),  flmn  =  f  (&m„)  ,  =  (&,£m,£n)  €  D. 

l,m,n— 0 

A  typical  collocation  approach  [16]  in  semi-discretizing  (3.1)  is  as  follows.  Find 
q  €  Vd  such  that 


E 


D.AC 


iNe  liNAJe)xNAQe) 


dqe 

dt 


\>e  dr 


V  r  ■  In  e 


ID,N. 


(»"') 


pedr 


n  ■ 


/  dD,Ne 


1 ; 


■Ne 


((3qe)*)  -INe  (?qe)]  • Pedr 


=  V  f  INc  (je)  lNe  (0e))  •  Pe  dr,  Vp  G  v„, 
„  J  D,Ne 


(3.2) 


where  Xjve  j  =  I nc  (In€  ( Jeal )  •  Inc  (30)  ■  The  direct  consequence  of  the  above 
collocation  is  that  the  integrand  in  each  integral  is  at  most  of  order  2Ne  in  each 
direction  ri,i  =  1,  2,  3.  The  subscript  Ne  in  the  integrals  means  that  the  integrals  are 
numerically  evaluated  using  the  corresponding  A^th-degree  LGL  (or  LG)  quadrature 
rule. 

4.  Semi-discrete  stability  for  conforming  meshes.  In  this  section,  we  pro¬ 
vide  a  stability  proof  for  the  both  elastic-acoustic  and  electromagnetic  cases  on  con¬ 
forming  meshes.  By  conforming  meshes  we  mean  that  the  intersection  of  two  elements 
is  either  an  entire  face,  and  entire  edge,  or  a  corner,  and  that  the  solution  order  is  the 
same  for  all  elements.  It  is  sufficient  to  prove  semi-discrete  stability  since,  by  a  result 
in  [19]  (and  the  references  therein),  if  the  semi-discrete  equation  is  stable,  the  fully 
discrete  system  with  the  time  derivative  discretized  by  a  locally  stable  Runge-Kutta 
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(for  example  the  classical  4th-order  Runge-Kutta  method)  is  stable  as  well,  as  long 
as  the  time  step  is  sufficiently  small. 

Here,  we  employ  the  energy  approach  to  prove  stability.  For  the  elastic-acoustic 
case,  the  semi-discrete  energy  functional  £d{t)  is  defined  as 

1  r 

£d(t)  ■=  E  £%e(t)  where  £eNe(t)  :=  -  /  (E  :  (C E)  +  pv  ■  v)  dx, 

e=l  Z  JO ‘,Ne 

and  for  the  electromagnetic  case  as 

-N el  -j  « 

£d(t)  :=  V  £%,  ( t )  where  £%.  (i)  :=  -  /  (eE  ■  E  +  p iH  ■  H)  dx. 

e=r  2  Jdb’n° 

For  convenience,  we  define  the  element-wise  discrete  L 2  inner  product  and  the 
induced  norm  on  a  generic  domain  D ,  which  could  be  an  element  or  its  boundary,  as 


(•bP  )d,n=  q-P  dx,  \\qfDN=  q-q  dx, 

Jd,n  Jd,n 

and  their  continuous  counterparts  as 

(q- P)d=  q -pdx,  \\qfD=  q-p  dx. 

J  D  J  D 

The  discrete  global  L2  norm  is  computed  as  the  summation  of  the  element-wise  con¬ 
tributions 

ni^  =  Eni2d^- 

e 

theorem  1  (Stability  for  conforming  meshes).  Assume  the  mesh  is  affine  and 
conforming  with  solution  order  N,  then  the  DGSEM  discretization  is  stable  in  the 
following  sense: 


2 

T>Ndd 


Moreover,  if  g  =  0,  then  j-t£d  <  0. 


Proof.  Substituting  p  := 


CE 

v 


for  the  elastic-acoustic  coupling  case, 


and  p  =  q  for  the  electromagnetic  case,  into  (3.2),  and  using  a  discrete  integration  by 
parts  formula  [24,  17],  we  obtain 


'D.N  L 


IN 


Vr  •  pe  —  Vr  •  In 


(3q! 


/  dD,N 


pe  dr 


dr 

JeQe  ■  pe  dr,  (4.1) 


'  D.N 


where  we  have  dropped  the  interpolation  operator  I nc  in  the  last  two  terms  on  the 
right  side  of  (4.1)  since  the  interpolation  is  an  orthogonal  projection  with  the  discrete 
L2  inner  product  [7],  e.g., 


/  D.N 


IN  ( q)pdx . 


(4.2) 


T.  Bui-Thanh  and  O.  Ghattas 


For  affine  meshes,  the  metric  terms  Jeal,i  =  1,2,3  are  constant,  and  thus 

Tn(t)  =Jeai-lN(d). 

After  some  simple  manipulations,  for  either  the  elastic-acoustic  case  or  the  electro¬ 
magnetic  case,  one  has 

J  (?n  (3qe)  •  Vrpe  -  Vr  ■  1N  (3qe)  •  pe)  dr  =  0. 

As  a  result,  the  volume  terms  vanish  and  the  preceding  equation  can  be  expressed  in 
the  integrals  over  physical  space  as 


E  /  - 

p  J  dD  e,N 


(5qe)*  -  x3qe 


pe  dx  +  /  ge 

J  D‘,N 


pe  dx.  (4.3) 


Next,  if  <9De  (~l9De  is  a  non-empty  two-dimensional  intersection  (<9De  fl3De  has  non¬ 
zero  two-dimensional  Lebesgue  measure)  we  combine  the  integrands  of  the  surface 
integrals  on  both  ”  and  “+”  LGL  (or  LG)  points.  This  is  possible  due  to  the 
mesh  conformity,  i.e. ,  the  number  of  LGL  (or  LG)  points  on  ”  and  “+”  sides  are 
the  same  and  they  can  be  reordered  to  be  exactly  coincident.  After  some  algebraic 
manipulations  for  the  surface  integrals  on  the  right  side  of  (4.3),  the  following  holds 
for  the  elastic-acoustic  case: 


X!  jt£d  =  -\  E  Jgo ,  N k°  {(n  ■  i61!)2 + p  cp p+ct  w2} 

+  fci  |||n  x  (n  x  [5])||2  +  +p~c~p+cf  || n  x  (n  x  [n])||2j  dx 


where  terms  involving  k±  are  zero  for  a  face  either  on  or  adjacent 
Similarly,  for  the  electromagnetic  case 


+  f  0e  •  P e  dx, 
to  the  acoustic  side. 


e!«=e /  v‘  v‘d* 

e  M  e  d  De  ,AT 

„  II"  x  "  x  Mil2  +  2pj) »"  x  "  x  lHIH2  dx’ 

where  terms  involving  E  and  H  vanish  for  PMC  and  PEC  boundaries,  respectively, 
and  both  vanish  for  dielectric  boundaries. 

Now,  for  0  =  0  the  overall  energy  is  non-increasing,  i.e., 


For  0^0,  we  obtain  the  following  estimate,  by  Cauchy- Schwarz  and  Young  inequali¬ 
ties, 


0  '  P  < 


2  (£d  +  \\^Ng\\ 


□ 
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5.  Mortar-based  non-conforming  approximations. 

In  this  section,  we  employ  the  mortar-based  approximation  idea  proposed  by 
Kopriva  et  al.  [14,  18].  However,  we  propose  to  use  a  discretized  mortar-based  ap¬ 
proximation  to  show  stability  in  addition  to  the  outflow  condition  and  global  con¬ 
servation  as  required  by  the  original  mortar  method  [14,  18].  As  will  be  shown,  our 
discrete  version  requires  a  special  quadrature  rule  in  order  to  simultaneously  satisfy 
all  the  requirements.  We  provide  a  setting  that  allows  unified  proofs  that  are  valid 
for  both  functional  (due  to  order  refinement)  and  geometric  (due  to  local  subdivision) 
non-conforming  approximations. 

The  following  conventions  are  adopted.  We  use  bold  face  type  to  denote  vectors 
of  nodal  values  of  the  corresponding  quantities  under  consideration.  For  example,  q 
is  the  vector  of  nodal  values  of  q.  In  addition,  we  use  upper  case  script  type  to  denote 
matrices,  e.g.  & . 

We  consider  non-conforming  approximations  due  to  domain  subdivision  in  which 
elements  may  be  subdivided  locally  while  their  neighbors  may  not.  For  simplicity 
of  exposition,  we  further  restrict  ourselves  to  the  case  where  a  subdomain  interface 
between  two  adjacent  elements  (two  elements  are  said  to  be  adjacent  to  each  other  if 
their  boundary  intersection  has  non-zero  two-dimensional  Lebesgue  measure)  must  be 
an  entire  face  of  either  of  them.  Nevertheless,  adjacent  elements  are  allowed  to  have 
different  solution  orders,  and  hence  order  refinement  (i.e. ,  p-refinement)  in  addition 
to  domain  subdivision  (i.e.,  /i-refmement)  is  permissible.  From  here  on,  by  “non- 
conforming  interface”  we  mean  that  an  entire  face  of  one  element  is  also  a  union  of 
faces  of  other  adjacent  elements  (/i-non-conforming),  or,  the  solution  orders  of  two 
elements  sharing  a  face  are  different  (p- non-conforming) .  A  non-conforming  interface 
with  one  element  on  one  side  and  seven  elements  on  the  other  side  is  shown  in  Figure 
5.1. 


Fig.  5.1.  A  non- conforming  interface  with  one  hexahedron  on  one  side  and  seven  hexahedra 
on  the  other  side. 

Consider  a  non-conforming  interface  where  on  the  ”  side  is  face  feo  of  ele¬ 
ment  e0  and  on  the  “+”  side  are  faces  fe.  of  elements  ei:i  =  l,...,Na.  Clearly, 
this  setting  includes  both  kinds  of  non-conforming  interfaces.  We  create  Na  mortars 
Mi,  i  =  1 ,Na  whose  ”  sides  are  seen  by  element  eo  and  “+”  sides  by  elements 
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e,;,  respectively.  As  in  [18],  the  geometric  order  on  the  mortars  must  be  the  lowest  geo¬ 
metric  order  of  the  contributing  elements,  and  the  polynomial  should  be  defined  along 
face  feo .  This  will  ensure  that  the  mortars  match  sub-interfaces  between  elements  ex¬ 
actly,  and  hence  the  metrics  Jeiaei  on  a  mortar  and  the  corresponding  contributing 
element  faces  are  represented  by  identical  polynomials.  The  solution  orders  of  the 
mortars  are  chosen  as  Nmi  >  max{7Veo, Nei]  which  is  sufficient  to  satisfy  the  outflow 
condition,  as  we  shall  show.  An  example  with  seven  mortars  corresponding  to  the 
non-conforming  interface  in  Figure  5.1  is  shown  in  Figure  5.2. 


Fig.  5.2.  Seven  mortars  corresponding  to  the  non- conforming  interface  in  Figure  5.1. 

The  goal  of  the  mortar  approximation  is  to  compute  the  contravariant  fluxes  on 
faces  fei ,  i  =  0, . . . ,  Na.  This  is  a  three-step  process.  First,  the  states  on  eo  and  e,  are 
projected  on  the  mortars  through  L2  projectors  and  i  =  1, . . . ,  Na- 

The  projected  states  are  then  used  to  compute  the  mortar  contravariant  fluxes  as  if 
each  mortar  is  a  conforming  face.  The  final  step  is  to  project  the  contravariant  fluxes 
back  on  the  element  faces  using  projectors  an(j  The  illustration 

of  steps  1  and  3  can  be  seen  in  Figure  5.2.  The  components  of  each  step  are  now 
detailed. 

State  qmi  on  the  ”  side  of  mortar  Mi  is  the  least  squares  projection  of  state 
qe°  from  element  eg  onto  space  7hvm. ,  be., 


WT(r)eVNm.  (5.1) 


Since  our  main  goal  is  to  prove  semi-stability,  least  squares  projection  of  the  type 
(5.1)  is  approximated  using  quadratures.  Nevertheless,  we  do  not  wish  to  violate  the 
outflow  condition  and  global  conservation.  In  fact,  the  outflow  condition  is  necessary 
to  ensure  stability,  as  shown  in  §6.  In  this  paper,  the  following  quadrature  rule  is  used. 
For  any  surface  integral  in  the  least  squares  projection  of  the  type  (5.1),  the  quadrature 
rule  is  chosen  to  correspond  to  the  integrand  with  the  highest  order.  For  example, 
the  lVmith-order  quadrature  rule  is  used  for  both  integrals  in  (5.1).  Explicitly,  we 
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approximate  (5.1)  as 


qm*7  (r)C*  M  dr  = 

qe0  ((Xe0)-l  Q  (p))  ^  (p)  ^  ygn*  (p)  g  ,  (5.2) 

which  yields 

qm«“  =  (^rmi)_1^e°^mi  qe° 
where  the  matrices  jfa /mi  and  f%e° _>mi  are  defined  as 

^2=  /  CWCW*. 

=  f  C  (r)C  ((^T1  o xe*  (r))  dr ,  wq  g  pN 

Similarly,  the  least  squares  projection  of  state  qe'  from  element  e%  onto  Vm.,,,  is  state 
qrrU  on  the  “+”  side  of  mortar  M. , .  Using  the  above  quadrature  rule  we  have 

q™+  (r)  C*  (r)  dr  =  [  qe-  (r)  (r)  dr,  W™'  W  G  ^ ,  (5.3) 

or,  equivalently, 

q™t  =  (jr1*)'1  S'*-*”1*  cf*  (5.4) 

V - V« - / 

with  defined  as 

=  [  IT  (0  q  (r)  dr,  Mq  G  VNf, . . 

Since  the  fluxes  depend  on  A,  /z,  c,  and  /.t,  we  perform  the  above  least  squares  projection 
procedure  on  them  as  well.  Based  on  the  projected  states  qm<  ,  qm<  ,  and  projected 
coefficients,  we  compute  the  contravariant  Riemann  fluxes  3*rl.  =  (j1 '  3(qm<))  on 
the  mortars  as  if  the  mortars  are  conforming  faces.  This  is  done  using  a  simple  relation 
between  the  contravariant  and  covariant  Riemann  fluxes  as  in  [16].  The  projected 

states  are  also  used  to  compute  the  contravariant  fluxes  3“ .  =  n-3  (qm;  ) ,  and  3m;  = 
h  ■  3  j  •  The  final  step  is  to  project  the  mortar  contravariant  fluxes  3*(. ,  3m„;  and 
3+.  onto  Vnb.  ,  i  =  0, . . .  ,Na.  Since  the  procedure  is  the  same  for  3*,, ,  3™  and  3+ 
we  describe  only  the  process  of  projecting  3m,.  to  obtain  the  contravariant  fluxes  3e; 
on  faces  /ei  of  the  contributing  elements.  The  discretized  least  squares  projection 
using  the  quadrature  rule  discussed  above  for  face  feo  is  as  follows: 


/  MijNrr 


[  3e0  (r)£f(r)dr  = 

d  feo  ’Ne0 

Na  r 

E  JM  N  K  (r)  q  ((X60)-1  o  xe-  (r))  dr,  Vq  G  VNeo ,  (5.5) 
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for  which,  in  matrix  notation,  the  vector  of  nodal  values  of  S'*  ,  F*  ,  is  computed  as 


Na 


^eoK  =  E- 


pnii—teof?* 


i= 1 


or  in  terms  of  projection  matrices, 

Na 


K  =  E^eo)-1' 


?>mi— >eo  rp* 


where 


=  [  £1°  (r)£j°  (r)  dr, 


>Mi,Neo 

Similarly,  the  vector  of  nodal  values  of  3t0  >  F“  giyen  by 

Na 


p-  _  \  '  Cprrii^e op- 
eo  /  _j 


(5.6) 


(5.7) 


(5.8) 


The  projection  to  compute  contravariant  fluxes  3)7  on  surface  fei  i  =  1  ,...,Na  of 
other  contributing  elements  is  simpler, 


K  (r)t?  ( r )  dr  = 


J  Mi, Ne. 

which  yields 


JMi^rr 


3*  (r)*?(r)dr  =  0,  Ve(r)e\,  (5.9) 


F*  = 


imij-te;  tti* 


where 


Similarly, 


=[  £V  (r)  £f  (r)  dr,  ^ro^ei  = 


Ft  = 


(5.10) 


(5.11) 


(5.12) 


Recall  that  the  outflow  condition  means  the  invariance  of  a  polynomial  function 
when  projected  to  the  mortars  and  back  to  the  face  [14,  18].  We  are  now  in  a  position 
to  discuss  the  outflow  condition  for  the  above  discrete  mortar-based  approximation. 

PROPOSITION  1  (Outflow condition) .  Assume  Nmi  >  max { Neo ,  JVeJ,  i  =  1 ,Na. 
Then  the  discrete  mortar-based  approximation  with  LG  quadrature  satisfies  the  strong 
outflow  conditions,  namely, 


rrii—tei  (fl0ei^rni  _ 


Opmi^ei  op 


=  I,  i  =  1 , . . . ,  Na , 


(5.13) 


Na 


o  ope^rrii  _  j 


(5.14) 
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where  I  is  the  identity  matrix  of  appropriate  size.  On  the  other  hand ,  discretizations 
using  LGL  quadrature  satisfy  the  outflow  condition  in  the  following  weak  sense, 


qei  ( r)£dr  = 


q Bi(r)£dr,  W  G  VNe 


(5.15) 


and 


J  qe°  ((^e°)_1  o  Xe>  (r))  £  ((Xe°)_1  o  Xe'  (r))  dr  = 


qe°  (r)  t  (r)  dr, 


WeVNeo,  (5.16) 


where  qEi  is  the  result  from  the  projection  of  qei  to  A ii  and  back  on  fei,  and  qe°  the 
result  from  the  projection  of  qe°  to  A 1,;  and  back  on  feo. 

Proof.  We  first  show  (5.13).  Denote  qm»  as  the  projection  of  qe;  on  mortar  A 1;, 
using  (5.3)  and  (5.9),  we  have 


I  Mi  ,Nn 


qmi  (r)  £  (r)  dr 


/Mi,N„ 


q'  (r)  £  (r)  dr 


/Mi, Nr, 


qBi  (r)  £  (r)  dr, 


lMi,Ne 


c\Bi  (r)  £  (r)  dr, 


W  G  VNmi , 

V£  G  Vnc  .  C  VNm .  ? 


which  imply  the  weak  outflow  condition  (5.15).  This  weak  outflow  condition  is  valid 
for  both  LGL  and  LG  quadrature  rules.  For  LG  quadrature,  however,  the  weak 
outflow  condition  also  implies  the  strong  one,  namely,  qe;  =  qe%  and  hence  (5.13), 
since  LG  quadrature  is  exact. 

For  (5.14),  denote  q as  the  projection  of  qe°  onto  mortars  Ai;.  Equations  (5.2) 
and  (5.5)  imply,  for  *  =  1, ... ,  Na, 


/  Mi,N„ 


qmi  (r)  £  (r)  dr  = 


f  qe°  ((Xe°)_1  o  XCi  (r)')  £  (r)  dr,  WeVNm., 


qmi  (r)£(  (X60)"1  o  (r) )  dr  = 


J  Mi, N, 

Na 

^  J  Mj  ,Nm 

[  qe°  (r)  £  (r)  dr,  V£  G  VNeo  ■ 
•/fen  We0 


(5.17) 


(5.18) 


Since  Nm.  >  ma x{Neo,  Ne.},  and  hence  Vnbo  C  ,  i  =  1, . . . ,  Na,  we  take  £  (r)  = 
£  ((r^'o!6’  (r))  in  (5.17)  and  sum  over  i  =  1, . . . ,  Na,  and  finally  subtract  from 

(5.18)  to  obtain  the  weak  outflow  condition  (5.16).  Again,  by  the  exactness  of  LG 
quadrature  we  have  qe°  =  qe°,  and  hence  (5.14).  □ 

remark  1.  It  is  clear  that  if  Nm.  =  Neo  =  Ne.,  i.e.,  only  h  non-conformity  is 
considered,  the  weak  outflow  conditions  are  indeed  strong. 
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PROPOSITION  2  (Global  conservation).  The  discrete  mortar-based  approximation 
satisfies  the  following  global  conservation, 


Na 

ca  =  E 

2=1  * 


K  dr. 


Proof.  This  result  is  an  easy  consequence  of  (5.5)  and  (5.9).  Indeed,  taking 
te°  =  1  in  (5.5)  and  =  1  in  (5.9),  we  have 


and 


Na 


(r)  dr  =  J2 

i=1  JMiiNrr 


Na 


Na 


E 


Mi, Nc 


K  ( r )  dr  =  E 


3*m,  (r)  dr, 


Mi,N„ 


5m,  {r)  dr. 


Observe  that  (J*o  (r)  is  a  polynomial  of  order  at  most  Neo,  while  £?*.  ( r )  is  a  polynomial 
of  order  at  most  Ne..  The  exactness  of  either  LGL  or  LG  quadrature  completes  the 
proof.  □ 

6.  Semi-discrete  stability  for  non-conforming  meshes.  We  start  this  sec¬ 
tion  with  a  discussion  on  why  the  outflow  condition  is  necessary  for  the  stability  proof 
to  hold.  The  semi-discrete  form  (3.2)  is  also  applied  for  non-conforming  approximation 

using  the  mortar  method  in  §5.  However,  the  contravariant  fluxes  ^5qe^  and  3qe  on 
the  boundary  <9D  are  the  projected  values  of  the  mortar  contravariant  fluxes  instead 

.  .  Accordingly,  using  commutativity  of  quadrature  and 
,  (3.2)  becomes 


of  the  trace  of  the  flux,  5qe 
integration  by  parts  [24,  17 


E 


'D,Ne 


dae 

jeQeJL.pedr- 


n 


'dD,Ne 


dt 

M 


—  5qe 


'  D,Ne 

-  5qc 


’-N, 


OD 


(5q' 


Vr  •  pe  dr 


P6  dr  =  E 


Jege  •  pe  dr. 


D,NP 


(6.1) 


conforming  approximation,  one  has  Sq11  .  ~  5qe  =  0,Vr  £  dD,  and  this  is  the 

1 9D 

reason  why  (4.1)  holds.  This  no  longer  holds  for  the  non-conforming  approximation 


unless  5qe  =  j?qc 


dD 


,  Vr  £  dD,  which  is  true  if  the  outflow  condition  is  satisfied 


and  Jqe  £  Vn A  sufficient  condition  for  (5qe  £  Vn‘  to  hold  is  that  the  medium 
properties,  i.e. ,  A,  p,  and  e,  are  element-wise  constant. 

We  introduce  the  global  interpolation  operator  n  whose  restriction  on  element 
De  is 


H  |  De  =  iNb, 

which  allows  us  to  obtain  the  following  stability  proof  for  our  non-conforming  approx¬ 
imation. 

theorem  2  (Stability  for  non-conforming  meshes  with  LG  quadratures).  As¬ 


sume 
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(i)  The  discrete  mortar  approach  in  §5  is  used  for  non- conforming  approxima¬ 
tions. 

(ii)  The  mesh  is  affine. 

(iii)  The  LG  quadrature  is  used,  i.e.,  the  strong  outflow  condition  is  satisfied. 

(iv)  A ,  /x,  and  e,  are  element-wise  constant. 

Then  the  DGSEM  discretization  is  stable  in  the  sense  that 

d  <  2  (Zd.  +  lln0llp«e1>d)  ■ 

Moreover,  if  q  =  (X  then  <  0. 

Proof.  As  discussed  above,  assumptions  (iii)  and  (iv)  imply  the  identity  5qe  „  — 

db 

5qe  =  0 ,  Vr  £  dD.  Therefore,  similar  to  the  proof  of  Theorem  1,  we  substitute 
for  the  elastic-acoustic  case  and  p  =  q  for  the  electromagnetic 

case,  to  obtain 


•  p edr+  [  Jege  •  pe  dr. 

JD  ,Ne 


(6.2) 


We  divide  the  surface  integrals  into  two  groups,  namely,  surface  integrals  associated 
with  conforming  and  with  non-conforming  interfaces.  The  former  group  has  been 
shown  to  be  non-positive  as  in  the  proof  of  Theorem  1.  We  therefore  need  to  con¬ 
sider  only  the  latter  group  for  which  we  take  a  typical  non-conforming  interface  and 
its  contributing  surface  integrals  as  in  §5.  We  shall  show  that  the  non-conforming 
contribution  is  also  non-positive. 

The  surface  integral  contributed  from  element  eo,  with  the  contravariant  fluxes 
projected  from  mortars,  can  be  written  as 


f/e0  ,JVeo  [£*e0  -  Ro]  • Pe°dr 

=  —  (Pe°)T  Me° 

using  quadrature 

=  -E&1  (P eof^m^eo  [f^  -  if-; 

using  (5.6)  and  (5.8) 

= -  (pm*_)T  [fe  - 

using  (5.3)  and  (5.7) 

=  -  £tai  [K  -  IK]  ■  pmr  dr 

using  quadrature 

=  -ZZ\fMi,Nmin-  [(3(qm*))*-|S  (qm,~  )  j  •  Pm’_ 

dr.  by  definition 

For  each  contributing  element  ei,i  =  1, . . . ,  Na,  the  surface  integral  on  face  fei 
with  the  contravariant  fluxes  projected  from  mortars,  is 

3e,  -  h&i  -Pe,dr 

=  -(Pe‘)TMe'  [F*.  -  if-J 

using  quadrature 

=  -(pe')T^"™  [F^-iF+; 

using  (5.10)  and  (5.12) 

=-(P  [F£-if+; 

using  (5.4)  and  (5.11) 

=  [K~hK]  -K  dr 

using  quadrature 

=  [(w*)]  -5®  (qm,+)]  -pmt  dr. 

by  definition 
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Therefore,  we  have  shown  that  each  non-conforming  interface  consisting  of  faces 
f£i  of  contributing  elements  e*,  i  =  0, . . .  ,Na,  is  equivalent  to  2 Na  conforming  faces 
associated  with  Na  mortars  Mj,  j  =  1, . . . ,  Na.  That  is,  the  surface  integrals  in  (6.2) 
in  fact  consist  of  conforming  interfaces — either  the  original  conforming  interfaces  or 
equivalent  conforming  mortars.  The  stability  proof  of  Theorem  1  for  conforming  faces 
can  now  be  applied  to  complete  the  proof.  □ 

When  the  strong  outflow  condition  is  not  satisfied,  i.e.,  when  LGL  quadrature  is 
used  as  in  Proposition  1.  We  have  the  following  stability. 

THEOREM  3  (Stability  for  non-conforming  meshes  with  LGL  quadratures).  Sup¬ 
pose  all  assumptions  in  Theorem  2  hold  except  that  the  LGL  quadrature  is  employed. 
Then  the  DGSEM  discretization  is  stable  in  the  sense  that 

~^Sd  <  -  ((1  +  2c) Ed  +  ||nfl||pjveiidj  , 


where  the  small  constant  c  converges  to  zero  if  the  exact  solution  q  is  sufficiently 
smooth. 

Proof.  When  the  numerical  integration  is  not  exact,  we  have  the  following  extra 
term 


E 


n  • 


•  pe  dr , 


which  can  be  shown  to  be  small  as 


E 


j  *. 

~Me 

.  ~3qe' 

>db,Nc 

3D  J 

■  pe  dr  <  c£d- 


(6.3) 


(6.4) 


We  now  provide  an  explicit  estimate  for  the  constant  c  to  show  that  c  is  indeed 
negligible.  It  is  sufficient  to  estimate  the  error  for  contributing  element  whose  face 
f£i  is  on  the  “+”  side  of  a  mortar.  From  the  weak  outflow  identity  (5.15)  and  the 
triangle  inequality,  we  have 


fl  • 

Mei 

.  -  3qe; 

•  pe‘  dr 

,iVe. 

3D  J 

fl 

•3qe‘ 

.  •  pBi  dr  —  [  \ 

li,Nm 

1 3D 

J  Mi 

< 


3D 


■  pe‘  dr 


+ 

/  h  •  dqei 

■  pe'  dr  - 

f  fi  ■  3qei 

■  pe i  dr 

JMi,Ne. 

3D 

J Mi 

3D 

Note  that  both  terms  on  the  right  side  of  the  preceding  inequality  are  of  the  same 
type,  namely,  the  error  between  LGL  quadrature  integration  and  exact  integration 
for  polynomial  of  order  2 N£i .  Since  Nmi  >  N£i ,  we  need  to  estimate  only  the  second 
term.  Using  an  error  estimate  from  [1]  together  with  Stirling’s  formula  we  have 


f  n  •  dqei 

■  pe*  dr  - 

[  fi  ■  5qei 

■  pe’  dr 

JMi,Ne. 

3D 

J  Mi 

3D 

=  O 


Nei+  1  (Ne,  -1 


2  Ne.  +  1 


N£, 


AN.. -2 


1  ,  2 Ne. 

tie 


-)2N. 
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which  shows  that  the  error  introduced  by  having  the  weak  outflow  condition  (instead 
of  the  strong  one  for  which  the  additional  error  (6.3)  is  zero)  is  decaying  exponentially 

with  respect  to  the  solution  order  Nei,  and  at  the  rate  2 Nei,  i.e. ,  he  ' ,  with  respect 
to  he.  Compared  to  the  optimal  error  rate  with  respect  to  he  and  Nei  in  (7.1)— (7.2), 
the  error  rate  resulting  from  the  weak  outflow  condition  is  much  smaller.  □ 

remark  2.  By  virtue  of  both  h  and  p  estimates  above,  we  see  that  even  though 
the  inexactness  of  the  LGL  quadrature  violates  the  strong  outflow  condition  and  hence 
generates  additional  boundary  terms  when  integration  by  parts  is  performed,  these 
terms  are  negligible.  From  now  on,  if  LGL  quadrature  is  used,  we  implicitly  absorb 
the  additional  error  terms  arising  from  the  weak  outflow  condition  into  constant  in 
the  convergence  estimate,  as  in  Theorem  f,  or  simply  ignore  them. 

7.  Convergence  rate  analysis.  The  previous  sections  show  that  our  discrete 
approximations,  both  conforming  and  non-conforming,  are  stable.  Together  with  con¬ 
sistency,  to  be  shown  below,  our  approximations  are  convergent  by  the  Lax-Ritchmyer 
equivalence  theorem.  The  direct  consequence  is  that  the  solution  can  grow  at  most 
exponentially  in  time,  which  is  typically  for  Lax-Ritchmyer  type  of  convergence.  This 
kind  of  convergence  result  is  interesting  for  theoretical  analysis,  but  may  not  be  appro¬ 
priate  for  assessing  the  actual  convergence  rate  of  a  numerical  method.  Fortunately, 
Hesthaven  and  Warburton  [12]  show  that  a  direct  convergence  analysis  is  possible, 
allowing  a  much  better  upper  bound  on  the  error  estimate.  We  adapt  this  type  of 
direct  convergence  analysis  to  derive  a  priori  error  bounds  for  our  non-conforming 
approximations . 

Recall  that  interpolation  introduces  truncation  and  aliasing  errors  [7,  16],  and 
hence  interpolation  is  generally  different  from  projection,  which  has  only  truncation 
error.  For  sufficiently  smooth  functions,  however,  the  aliasing  error  either  is  spectrally 
small  [7,  11,  16]  or  can  be  made  equal  to  zero  [12].  Following  [12],  we  shall  make  no 
distinction  between  interpolation  and  projection  in  what  remains. 

The  following  conventions  are  used  in  this  section.  We  reserve  q  for  the  unknown 
exact  solution,  Xve  q  for  the  projection  of  q  on  Vnc  ,  and  q  jve  the  solution  of  the  discrete 
form  (3.2)  restricted  on  element  e.  In  addition,  C  denotes  a  generic  constant  that  may 
have  different  values  in  different  contexts,  q^  is  defined  by  q<j|De  =  <live)  and  finally 
a  dummy  variable  q  lives  in  different  spaces  for  different  inequalities.  To  begin,  we 
recall  the  following  fundamental  hp  approximation  error  bounds  [2,  3,  4], 


\Q  -^Neg\\ /Da\  <  C 


h° 


NS 


HSe  ( De ) 


De) ,  0  <  r  <  se 


hee  1^2  1 

k  —  |laD=  —  C  s  _  1/2  lklliT<»e(De)  5  Se  >  A 


(7.1) 

(7.2) 


with  he  =  diam(  De),  cre  =  min{iVe  +  l,se},  and  |Hl.ff»qDe)  denoting  the  usual  Sobolev 
norm. 

For  approximation  using  tensor  product  LGL  quadrature,  the  following  equiva¬ 
lence  of  the  discrete  and  continuous  norm,  an  extension  of  the  one  dimensional  version 
in  [11],  is  useful  in  passing  from  the  discrete  norm  to  the  continuous  one  and  vice  versa: 
Vg  G  Vn, 


klip  —  Ikll  d.jv  — 


„  i  \3/2 
2+» 


ID  ’ 


and 


lab  — 


I  dD.N  — 


2  + 


N 


kll 


d  D  ' 
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We  first  derive  the  convergence  rate  for  conforming  meshes.  Since  the  electro- 
magnetic,  acoustic,  elastic,  and  coupled  acoustic-elastic  wave  equations  are  similar, 
we  analyze  the  electromagnetic  case  and  leave  out  details  of  the  others.  Denote 
Tq  =  [Te,Th]  as  the  the  truncation  error  that  results  from  substituting  the  exact 
solution  q  in  the  discrete  equation  (3.2).  Using  the  fact  that  q  satisfies  the  Maxwell’s 
equations,  we  have 

{^k,^NeTE)De  ^  =  (ek,xNy  x  (H  —  iNeH))De  Ne 

+  (tk,  0  ffL  n  x  (Z+  [XNcH]  -  n  x  [Z^-E]))  ,  V4  &VNe 

\  ^  1 1 ^  if  /  dDe,Ne 

(tk,INeTH)De  =  (£k,XNV  x  (E  —  INcE))Da  N^ 

+  ( £k ’  o  (mi n  x  ( Y+  \xNeE]  +  n  x  [lNeH])\  ,  V4  G  ZVe 

V  z  ll1  JJ  / 

where,  to  the  end  of  this  section,  p  and  e  are  assumed  to  be  element-wise  constant, 
and  the  mesh  is  affine.  Note  that  since  the  truncation  errors  for  the  electric  and 
magnetic  equations  are  similar,  we  need  to  estimate  only  the  former  and  infer  the 
later.  Since  I ncTe  G  Vnc,  we  can  take  4  =  XnbTe  and  use  the  Cauchy-Schwarz 
inequality  together  with  the  equivalence  of  discrete  and  continuous  norms  to  obtain 

llJ^||2De>JVe  <  272  II iNy  X  {H-1NH) ||De  ll^.r-11,,. 

+  92  x  (Z+  pNeH]  -nx  [INE ]) 

Now  using  the  following  inverse  trace  inequality  [23],  \/q  G  Vnc, 

UW<C^\\q\\De  (7.3) 

lie 

yields 

\\inte\\d^  <  C\\lNy  x  (H-lNeH) ||De 

+  CTU-2  7^i(Z+[lNBHr}-[INET])  ,  (7.4) 

he  ZUZ//1  dD‘ 

where  we  have  introduced  the  tangent  components  of  E  and  H  as 
Et  =  n  x  (n  x  E) ,  HT  =  n  x  H. 

We  now  have  the  following  consistency  result. 

LEMMA  1  (Consistency).  Suppose  that  each  component  qf  G  HSe  (De)  ,se  >  3/2, 
for  i  =  1, . . . ,  d,  with  d  =  6  for  electromagnetic  case  and  d  =  12  for  elastic-acoustic 
case.  Furthermore,  assume  that  p  and  e  (X  and  p  for  elastic- acoustic  case)  are 
element-wise  constant,  and  the  mesh  is  affine.  There  exists  a  constant  C  dependent 
on  s,  angle  condition  of  De,  and  local  values  of  p  and  £  (X  and  p  for  elastic-acoustic 
case),  but  independent  ofq,he,  and  Ne  such  that 

\\XNsTq\\VNeltd  <  C  ^  _ 3/2  l!Clll[ffse(D'=)]d  • 

e  -Ne 
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Proof.  Since  the  proofs  for  electromagnetic  and  elastic-acoustic  cases  are  similar, 
we  provide  only  the  proof  for  the  former.  We  begin  by  estimating  the  bound  for  the 
first  term  on  the  right  side  of  (7.4).  Using  approximation  result  (7.1)  we  have 

||I,veVx  (H-lNeH)\\De  <||Vx  {H-lNH)\\De  <Cj^\\H\\[H.a(De)]s. 

To  estimate  the  second  term  we  use  the  regularity  of  the  exact  solution,  i.e. ,  the 
tangent  component  of  the  field  is  continuous,  and  the  triangle  inequality  to  bound 
\\[lNeHT]\\dD„  and  hence  similarly  for  ||  [JjveST]||aDe,  as 

||[2Jve-HT]||aDe  <  \\H~T  -lNH-T\\dDB  +  \\Ht  ~lNH+T\\dDe 

Two  terms  of  the  right  side  of  the  preceding  inequality  are  of  the  same  type,  and  there¬ 
fore  we  need  to  estimate  only  the  bound  for  the  first  term.  But  this  is  straightforward 
using  (7.2),  i.e., 


L^e-l/2 


,  CTe  — 1/2 


\H-  -  1N,  H 


00*  <C 


N' 


e- 1/2 


l^r||[ff..(D.)]3<C' 


n;_ 


e_ 1/2 


\H\ 


[H‘e(D*)]3  ■ 


Now  combining  the  above  estimates  for  both  terms  on  the  right  side  of  (7.4),  summing 
over  all  elements,  and  using  the  discrete  Holder  inequality  completes  the  proof.  □ 

remark  3.  Note  that  the  proof  of  Lemma  1  is  carried  out  for  conforming  meshes. 
However,  by  the  proof  of  Theorem  2,  the  surface  integrals  in  (6.2)  in  fact  consist  of 
conforming  interfaces — either  the  original  conforming  interfaces  or  equivalent  con¬ 
forming  mortars.  Hence,  the  pi'oof  for  non- conforming  meshes  is  almost  identical 
except  for  bounding  the  boundary  terms  which  are  now  defined  on  the  mortars  instead 
of  on  the  contributing  element  faces.  But  the  fields  on  the  mortars  are  the  L 2  orthog¬ 
onal  projections  of  those  on  the  contributing  element  faces,  which  implies  L 2  norms 
of  fields  on  mortars  to  be  at  most  those  on  contributing  element  faces.  This  shows 
that  the  proof  for  conforming  meshes  is  sufficient. 

We  now  state  the  first  convergence  result. 

THEOREM  4.  Assume  qe  €  [ Hs *  (De)]d,se  >  3/2  with  d  =  6  for  electromagnetic 
case  and  d  —  12  for  elastic-acoustic  case.  In  addition,  suppose  qc/(0)  =  nq(0).  Fur¬ 
thermore,  assume  that  p  and  e  (X  and  p  for  elastic- acoustic  case)  are  element-wise 
constant,  and  the  mesh  is  affine.  Then,  the  solution  q d  of  the  discrete  form  (3.2) 
converges  to  the  exact  solution  q,  i.e.,  there  exists  a  constant  C  that  depends  only 
on  the  angle  condition  of  De,  s,  and  the  material  constants  p  and  e  (X  and  p  for 
elastic-acoustic  case)  such  that 


|q  (t)  -  qd(t)|| 


V™el,d 


K 

Nf 


|q  (^:)ll[H»e(De)]<i  +  ^ 


hZ° 

Ni 


o/o  max 
*~3/2  [o,t] 


I  [Hse  (De)]d 
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Proof.  We  begin  the  proof  with  the  following  identities 

.  d  E  —  ENe 


(In^E  —  ENe,lN(lTE)De  N  =  [INbE  —  E 


'■‘N.1  £- 


dt 


D=,JVe 


—  (ENeE  —  ENe,2NeV  x  (I NeH  —  HNJ)De  Nt 
1 


En€E  —  ENe,  2  {{z}}  n  X  ^ +  ^ Ne H  _  H N el  ^  n  X  PiVe-®'  —  ^jvJ) 


aDe,Afe 


(ljve  H  —  HNe,  lNe  T 


’D‘.NC 


—  (  ENcH  —  HNe,p 


d{lNH-HNe) 

dt 


De  ,Ne 


(IjVe-ff  —  HNs,INfS7  x  (IncE  —  ENJ)De 


,AT„ 


-  (lNBH-HNe, 


2{{Y}} 


n  x 


(y+  -  £jve]  +  n  x  [JWen  -  iTjvJ) 


9Dc,Are 


where  we  have  substituted  the  exact  solution  into  the  discrete  equation  (3.2),  and 
used  1ncE  —  Exe  and  InbH  —  H jve  as  test  functions  for  the  electric  and  magnetic 
equations,  respectively. 

Following  the  proof  of  Theorem  2,  we  integrate  the  preceding  equations  by  parts, 
sum  up  the  resulting  equations,  cancel  the  volume  terms  involving  the  fluxes,  and  use 
the  non-positiveness  of  the  boundary  integrals  to  arrive  at 

—  II Idq  —  C\d\\x>Nei,d  —  —  c|Are;2jve2nq)Dc,]Ve,  ) 

e 

where  we  have  used  the  fact  that  the  material  constants  /i  and  e  are  bounded  away 
from  zero.  Next,  we  use  Cauchy-Schwarz  and  then  the  discrete  Holder  inequalities, 
then  apply  the  consistency  result  of  Lemma  1  to  obtain 


^  linq  q<i |lx>Nei,d  <  c ^ 


-,(Te  —  1 


[ff=e(D=)]d  <  CY1  ATSC- 3/2  max||q||[/f8e(De)]ci  , 


N!b~3/2  . l"  "  ~  ^  [o,t] 


which,  after  integrating  in  time,  yields 


|nq  (t)  -  qd  (t)Hx,«eljd  <  Ct}2  N*e_ 3/2  max || q  (t)||[Hae(De)]<i , 


where  we  have  used  q<i(0)  =  nq(0).  Now,  using  the  triangle  inequality  we  have, 


llq  (t)  -  q d  it) \\VNel  d  <  |jq  (t)  -  nq  (t)||x,^el +  ||Hq  (t)  -  qd  (t)||x,«el,d  • 

Finally,  using  the  equivalence  of  the  discrete  and  continuous  norms  and  (7.1)  ends  the 
proof.  □ 

REMARK  4.  Since  all  norms  are  equivalent  in  finite  dimensional  spaces,  the  result 
of  Theorem  4  holds  for  other  norms  as  well,  with  possibly  different  constant  C .  In 
particular,  we  have 


53l|q(t)-qd(*)llD.lJVa<C53 


[HSe  (De )]‘ 


hz« 


-—7-  max 

e-3/2  [o.t] 


|q(*)ll 
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remark  5.  We  again  emphasize  that  the  LGL  quadrature  introduces  additional 
error  terms  when  the  discrete  integration  by  parts  is  performed.  However,  they  are 
shown  to  be  negligible  in  Theorem  3  and  Remark  2,  so  their  absence  in  the  proof  of 
Theorem  4  is  justified.  The  proof  for  the  LG  quadrature  follows  exactly  the  same  line 
and  hence  is  omitted  here.  The  only  difference  is  that  all  the  numerical  integrations 
are  exact,  and  therefore  there  is  no  need  to  invoke  the  equivalence  of  discrete  and 
continuous  norms. 

The  proof  of  Theorem  4  is  simple  since  it  directly  uses  the  consistency  and  stability 
results.  Nevertheless,  the  rate  is  suboptimal  in  h  by  a  factor  of  1/2  and  in  TV  by  1 
compared  to  the  DG  literature.  This  loss  is  incurred  when  estimating  the  truncation 
error  and  using  the  inverse  trace  inequality  (7.3).  This  begs  for  a  more  sophisticated 
proof  that  does  not  lead  to  a  deterioration  of  the  optimal  convergence  rate.  Inspired 
by  the  work  of  Warburton  [26] ,  we  present  a  proof  that  improves  the  convergence  rate 
in  h  by  a  factor  of  1/2  and  in  N  by  1.  We  start  by  rewriting  (3.2),  after  integrating 
by  parts,  for  an  affine  mesh  and  element-wise  constant  materials  as 

E  [  “Iw  •  P6  rf£C  ~  /  SqATe  •  (Va;  •  P6)  dx 
e  J  De  at  J  Dc 

+  [  n  •  [(SqjvJ*]  •  pe  dx  =  E  [  ge  ■  pe  dx +  Ke,  Vp  G  Vd,  (7.5) 
Jd  Dc  e  J  D« 

where,  by  the  proof  of  Theorem  3,  lZe  =  o  and  lZe  =  0  for  LGL  and  LG 

quadratures,  respectively.  Since  the  exact  solution  q  also  satisfies  (7.5)  with  lZe  =  0, 
we  obtain 


E 


d(q  -  qjve) 

dt 


P edx-  £(q  -  qjve)  •  (vx  •  pe)  dx 

J  De 

+  [  h  ■  (S'  (q  -  qjvJ)*  •  pe  dx  =  -  Y'  Tle,  Vp  e  Vd- 
Jo  D= 


(7.6) 


Next,  adding  and  subtracting  the  elemental  projection  of  the  exact  solution,  i.e.,  2/veq, 
yield, 


E 


d  (2/veq  -  qjve) 
dt 


pedx-  3(lweq  -  q^J  •  (Vx  •  pe)  dx 
J  De 

+  f  n-($  (INeq  -  q nJ)*  ■  pe  dx 
JdD‘ 

=  -E  f  ^{q-^N^q))* -pe dx-ne,  vpevd, 

e  J™* 

where  we  have  used  the  linearity  of  J  and  the  following  orthogonal  identity, 


(7.7) 


[  (q-lNeq)-pedx  =  0,  pe  £  L2  (De)d  . 

J  De 


Now,  integrating  (7.7)  by  parts,  testing  the  resulting  equation  and  (7.7)  with  p 
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Iiveq  —  qAre,  and  then  adding  them  yield, 


nq  qd|lx>«ei, 

d  +  £^/  «'(S(q-  INcq  ))*  •  [2)veq  -  qjve]  dx  +  Ke 

e  ZJ9°‘ 

=  -£/  «■ 

p  JdD‘ 

(S  (2/veq  -  qjve))*  -  -S (2/veq  -  qjvJ 

•  (InA  -  qjve)  dx 

_ s 

- v - 

Ke 

(7.8) 


On  the  other  hand,  if  p,  A,  p  and  e  are  bounded,  the  proof  of  Theorems  1  and  2  implies 

Ke>ae\\[lNeq-qNe]\\2[LHdD.r,  (7-9) 

for  some  ae  >  0. 

We  are  now  in  the  position  to  prove  the  following  convergence  result. 
theorem  5.  Assume  qe  €  [HSe  (De)]d,se  >  3/2  with  d  =  6  for  electromagnetic 
case  and  d  =  12  for  elastic-acoustic  case.  In  addition,  suppose  qd(0)  =  Ilq  (0) .  Fur¬ 
thermore,  assume  that  p  and  e  (X  and  p  for  elastic- acoustic  case)  are  element-wise 
constant,  and  the  mesh  is  affine.  Then,  the  solution  q^  of  the  discrete  form  (3.2) 
converges  to  the  exact  solution  q,  i.e.,  there  exists  a  constant  C  that  depends  only 
on  the  angle  condition  of  De,  s,  and  the  material  constants  p  and  e  (X  and  p  for 
elastic-acoustic  case)  such  that 


|q(*)-q<»(t)ll 


V”eld 


ftSe  II  q  (i)ll[H»e(D«‘)]d  +  * 


=  -1/2 


m 


— T7T  max 
‘-1?2  [o,t] 


l[ffse(De)]d 


Proof.  Using  the  Young  inequality  for  some  /3e  >  0  we  have, 

£  /  \™- ($(c\-Zncc\))* qNe}\  dx  <J2f3e\\q-lNeq\\2[L2{dDe)]d. 

+  £  llP^q  ~  q^jllfL^aD*)]^  >  (7-io) 

where  we  have  used  the  linearity  of  S*-  Combining  (7.8),  (7.9),  and  (7.10)  yields 
—  ||IIq  —  qd||i>Nei ,d  +  £  —  ll[-^veq  —  clA'e]ll[L2(c>D«)]d 

<  £/3e  Hq-^Ar.qllfi^aD*)]-  +  \Ke\  ■  (7-H) 

e 

Next,  taking  f3e  >  \  and  using  (7.2)  give 

^  _  p2ac  —  l 

-  ||nq  -  qd\\x)Nel  d  <CJ2  -^=1  max|jq||[ffSe(D.)]£i  • 

The  rest  of  the  proof  follows  similarly  to  that  of  Theorem  4.  □ 

We  have  restricted  ourselves  to  the  case  of  affine  hexahedral  meshes  and  element¬ 
wise  constant  medium  properties.  This  allows  us  to  eliminate  the  volume  integrals 
of  fluxes  on  the  right  side  of  equation  (4.1)  since  differentiation  and  interpolation 
commute.  The  discrete  stability  is  then  apparent  due  to  the  non-positiveness  of  the 
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surface  integrals  of  fluxes.  Most  of  the  results  still  hold  for  meshes  with  Ja‘  = 
constant,  for  i  =  1,2,3.  In  addition,  it  is  not  difficult  to  see  that  the  stability,  and 
hence  convergence,  results  for  both  conforming  and  non-conforming  approximations 
are  still  valid  for  meshes  with  curved  elements,  provided  that  the  contravariant  fluxes 
are  polynomials  with  order  at  most  Ne,  i.e. , 

S!  =  JV.5e^«,  *  =  1,2,3,  (7.12) 

for  which  we  again  have  the  commutativity  of  differentiation  and  interpolation.  To¬ 
gether  with  the  metric  identities  in  [15],  we  again  can  eliminate  the  volume  integrals 
after  integrating  by  parts.  Moreover,  it  is  clear  that  our  results  remain  true  for  other 
types  of  meshes,  e.g.,  affine  tetrahedral  meshes  as  well,  as  long  as  the  discrete  integra¬ 
tion  by  parts  is  possible  and  the  commutativity  of  differentiation  and  interpolation  is 
valid.  For  general  curvilinear  hexalredral  meshes,  it  is  not  clear  whether  the  volume 
integrals  vanish  (or  become  negative)  or  not  since  differentiation  and  interpolation 
generally  do  not  commute  (even  over- integration  is  not  helpful  in  this  case). 

8.  Conclusions.  We  have  presented  an  analysis  of  a  non-conforming  hp  discon¬ 
tinuous  Galerkin  spectral  element  method  for  time  domain  solution  of  wave  propaga¬ 
tion  problems  in  acoustic,  elastic,  coupled  elastic-acoustic,  and  electromagnetic  me¬ 
dia.  We  have  proven  consistency,  stability,  and  convergence  under  the  usual  assump¬ 
tions,  i.e.,  affine  meshes  and  element-wise  constant  medium  properties.  Our  analytical 
results  hold  for  both  exact  numerical  integration  using  tensor  product  Legendre-Gauss 
quadrature  and  inexact  numerical  integration  using  tensor  product  Legendre-Gauss- 
Lobatto  quadrature.  The  key  ingredient  of  our  proposed  approach  is  the  development 
of  a  discrete  mortar-based  approach  for  non-conforming  approximations.  With  this 
mortar  construction,  we  have  shown  that  the  proofs  for  non-conforming  cases  closely 
follow  those  for  conforming  ones,  and  the  resulting  convergence  retains  the  same  op¬ 
timal  rate  as  for  standard  DG  methods.  Numerical  experiments  in  [27]  for  acoustic, 
elastic,  and  coupled  acoustic-elastic  wave  propagation  on  h-non-conforming  meshes 
confirm  the  theoretical  rates  presented  here,  and  a  further  companion  paper  numeri¬ 
cally  demonstrates  optimal  convergence  for  electromagnetic  wave  propagation  [6]. 
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